#include "../include/Kernel.hpp"
#include "../include/Random.hpp"
#include <armadillo>
#include "../include/Density.hpp"
#include <iostream>
#include "../include/Distributions.hpp"
#include "../include/utils.hpp"

#define M 8


extern "C"{

void MHf(int *n, int *p, double *Y, double *pX, int *pM)
{
	mat X=vect2mat<double>(pX,n[0],p[0]);
// test Density
	Density::Density_Probit *DP= new Density::Density_Probit(Y,X);
	//cout << DP->Density(theta,3,4) << " ";	
	//cout << DP->Phi(1.96);
	//test Proposal
	mat m(p[0],1);
	mat foo(p[0],p[0]);
	m.fill(0);
	foo.fill(0.01);
	mat s=diagmat(foo);
//	cout << s;
	Distribution::Gaussian *G=new Distribution::Gaussian(p[0],m,s);
	mat theta=G->r();
	//for(int i=0;i<10;i++){cout << A(i,0) << " ";}
	cout << theta << "\\";
	Kernel_MH<Density::Density_Probit,Distribution::Gaussian> K(DP,G);
	mat sum=theta;
	for(int i=0;i<pM[0];i++){
	
		mat temp=theta.t();	
		K.Move(&(theta),n[0]-1);
		//theta=temp.t();
		sum=sum+theta;
		cout << sum.t()/i << "\n";	
	}

}

}
